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ABSTRACT 

We use hydrodynamical simulations from the OWLS project to investigate the dependence of 
the physical properties of galaxy populations at redshift 2 on metal-line cooling and feedback 
from star formation and active galactic nuclei (AGN). We find that if the sub-grid feedback 
from star formation is implemented kinetically, the feedback is only efficient if the initial 
wind velocity exceeds a critical value. This critical velocity increases with galaxy mass and 
also if metal-line cooling is included. This suggests that radiative losses quench the winds 
if their initial velocity is too low. If the feedback is efficient, then the star formation rate is 
inversely proportional to the amount of energy injected per unit stellar mass formed (which is 
proportional to the initial mass loading for a fixed wind velocity). This can be understood if 
the star formation is self-regulating, i.e. if the star formation rate (and thus the gas fraction) 
increase until the outflow rate balances the inflow rate. Feedback from AGN is efficient at high 
masses, while increasing the initial wind velocity with gas pressure or halo mass allows one 
to generate galaxy-wide outflows at all masses. Matching the observed galaxy mass function 
requires efficient feedback. In particular, the predicted faint-end slope is too steep unless we 
resort to highly mass loaded winds for low-mass objects. Such efficient feedback from low- 
mass galaxies (M» «: 10 10 M G ) also reduces the discrepancy with the observed specific star 
formation rates, which are higher than predicted unless the feedback transitions from highly 
efficient to inefficient just below the observed stellar mass range. 

Key words: cosmology: theory - galaxies: formation - galaxies: evolution - galaxies: fun- 
damental parameters - methods: numerical 



1 INTRODUCTION 

Simulating the growth of dark matter (DM) haloes from initially 
small density perturbations through to the present day has become 
well established. Even the complex, non-linear stage of structure 
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formation can be predicted by means of high-resolution gravita- 
tional ,/V-body simulations. The distribution of DM haloes derived 
from these simulations agrees very well with observations. The for- 
mation and evolution of galaxies is, however, much less well un- 
derstood. Modeling the baryonic component is much more difficult 
than simulating the DM due to the collisional nature of the gas and 
the wealth of phenomena that need to be taken into account (cool- 
ing, star formation, feedback, etc.). 

There are two popular approaches to tackle this challenging 
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task. In semi-analytic models, simple descriptions of the behaviour 
of the baryonic component, as a function of the DM halo mass, 
merging hist ory and environment, describe the evolution of gas 



and stars (e.g. Kauffmann et al.l 1 999 : So merville &Primack! l999: 



ICole et alE oOO). The freedom to choose functional forms and pa- 
rameter values combined with the ability to run large numbers 
of models, ensure that reproducing observations is usually within 
reach. While this approach has great advantages, such as the abil- 
ity to make mock galaxy surveys that are sufficiently realistic to 
reveal observational biases, there are also significant drawbacks. 
The large number of parameters can make it difficult to identify the 
key physical processes. More importantly, the ability to reproduce 
observations with a model that uses unphysical functional forms 
or unrealistic parameter values to describe physical processes can 
easily result in erroneous conclusions and misplaced confidence. 

The other approach is to follow both the dark matter (DM) 
and the baryonic components by direct simulation. While the DM 
is nearly always simulated using particles, baryons can either be 
modeled with Euleri an methods (dis c retizing the volu me in an 
(adaptive) grid, e.g. iRvu et alj 11999 ; ICen et alj Il990h or using 
the Lagrangian approach also used for t he DM (discretizing the 
mass using particles e.g. lEvrardl fl988l ; iHernquist & Katj [l989; 
Th omas & Couchmar]|l992h . Here, the freedom is limited to the 
parametrization of unresolved sub-grid processes. The high com- 
putational expense associated with full numerical simulations and 
the reduced level of freedom in the sub-grid modelling mean that, 
thus far, numerical simulations have been much less successful in 
reproducing observations of galaxy populations than semi-analytic 
models. Compared with the semi-analytic method, the advantages 
of the simulation approach include the much reduced (though still 
present) risk of getting the right answers for the wrong reasons, the 
ability to ask more detailed questions due to the tremendous in- 
crease in resolution, and the fact that not only galaxies, but also the 
intergalactic medium (IGM) is modeled. 

As many processes related to the baryons are not (well) re- 
solved by even the highest resolution simulations, they are dealt 
with in the so-called sub-grid models. Among these are radiative 
cooling, the temperature and pressure of the multiphase gas at high 
densities (in the rest of the paper loosely called "the interstellar 
medium (ISM)") and the formation of stars, the energy and mo- 
mentum fed back by these stars into the ISM/intra-cluster medium 
(ICM)/IGM, stellar mass loss and the growth of supermassive black 
holes and associated feedback processes. 

In this work, we employ cosmological, hydrodynamical sim- 
ulations to investigate a number of basic baryonic properties of 
haloes, including: the (specific) star formation rate (SFR); the stel- 
lar, gas and baryon fractions; the gas-consumption timescale and 
the galaxy stellar mass function. Reproducing observations in de- 
tail is not the main goal of this paper and so we have not attempted 
to tune our models or to optimize the sub-grid implementations to 
match any particular data-set. Instead, we focus on understanding 
how different physical mechanisms shape the galaxy population. 

We make use of a subset of the simulat ions from the Over - 
Whelmingly Large Simulations project OWLS dSchave et alfcoioh . 
a large (~ 50 simulations) set of cosmological, smoothed particle 
hydrodynamics (SPH) simulations, each of the same volume of the 
universe, run with a wide variety of different prescriptions for sub- 
grid physical processes. The large variety of input physics in the 
OWLS runs allows us to investigate properties of haloes and their 
relation to the physical and numerical parameters. In this paper, we 
investigate the effect of feedback from star formation and from ac- 
creting supermassive black holes, as well as metal-line cooling. In 



a companion paper faaasetalj|2013l , hereafter Paper II) we inves- 
tigate the effects of the assumed cosmology, the reionization his- 
tory, the treatment of the unresolved, multiphase ISM, the assumed 
star-formation law and the stellar initial mass function. These in- 
gredients are all present in the simulations presented here, but the 
parameters are not varied. 

This work ( together with Paper II) complements that of 
ISchaveetailfcOlOl) . where the cosmic star formation histories pre- 
dicted by the OWLS models were analyzed. The global SFR density 
can be decomposed into a DM halo mass function, which is deter- 
mined by the cosmology, and the statistical distribution of the SFR 
as a function of halo mass. Here we will study the latter, which is 
astrophysically more relevant than the global SFR density as it re- 
moves the main effect of cosmology (the mass function) and allows 
us to investigate how the effects of the various baryonic processes 
vary with mass. W hilst we will add a dimension to the work of 
ISchave et all d2010h by investigating the dependence on mass, we 
will remove another one in order to keep the scope of the study 
manageable. Thus, we limit ourselves to z = 2 and to the high- 
resolution series presented in ISchave et al.U20ld) (these runs were 
halted at this redshift). 

This paper is structured as follows: In Sec. [2] we describe the 
reference simulation used in this study (Sec. 12. it , how we define 
and extract haloes (Sec. 12.2b and how we compare the simula- 
tions to observations (Sec. 12.3b . In Sec. [3] we describe how the 
properties of the galaxies in our reference simulation depend on 
both halo mass (Sec. 13.1b and galaxy stellar mass (Sec. 13.2b . In 
Sec [4] we describe how galaxy properties depend upon the physics 
included in the simulation. In this paper we consider the effects 
of: metal-line cooling (Sec. 14. lb , constant-energy supernova (SN) 
winds (Sec. 14. 2t , the effect of a top-heavy IMF in high-density gas 
(Sec. 14. 3b , "momentum-driven" winds (Sec. |4.4b , winds decoupled 
from the hydrodynamics (Sec. 14.5b . thermally coupled SN-driven 
winds (Sec. 14. 6b . and the effect of feedback from active galactic nu- 
clei (AGN) (Sec. 14. 7b . Finally, in Sec.[5]we summarize our findings 
and conclude. In Appendix lAl we study the numerical convergence 
of our results using simulations of different volumes and with dif- 
ferent numerical resolution while in Appendix [B] we show that our 
results are insensitive to our particular choice of halo finder. 



2 NUMERICAL TECHNIQUES 

For a detaile d discussion of the f ull set of OWLS models we refer 
the reader to ISchave et alj ( bold) . Here we briefly summarize the 
reference model. Throughout this paper we refer to this reference 
simulation as 'REF\ 



2.1 The reference simulation 

We ran our simulations in periodic boxes of 25 co-moving ft^'Mpc 
with 512 3 dark matter and baryonic particles (which originally are 
collisional 'gas' particles, but can be converted into collisionless 
'star' particles in the course of the simulation). We evolved the 
particles using an extended version of the N-Body Tree/SPH code 
Gadget3 (last described in ISpringell2005h . The simulation particles 
have masses of 8.68 x 10 6 /r' M Q for dark matter and 1.85 x 10 6 /r' 
M G for baryons (initially; the baryonic particle masses change in 
the course of the simulation due to mass transfer from stars back 
to the gaseous phase). In Appendix lAl we show that our results are 
reasonably well converged with respect to the resolution and box 
size of the simulation. The gravitational softening length is initially 
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fixed at 1/25 the inter-particle spacing in co-moving coordinates, 
but below z - 2.91 it is fixed at 0.5 / r'kpc in proper units. 

CMBFAST (Seljak & Zaldarriaga 19961) was used to generate 
initial cond i tions, that were evolved forward in time using the 
IZerDovichl dl97d) approximation from an initial glass-like state. 
The simulation is started at z = 127. The value of the cosmolog- 
ical parameters are n m = 0.238, Si A = 0.762, £l b = 0.0418, h = 
0.73, erg = 0.74 and n s = 0.951. These valu es are consistent with 
the 7-yr WMAP data iKomatsu et alj|201lh . The variation in sub- 
grid physics implementations is the main power of the OWLS set of 
simulations. The rest of this paper deals with variations of the sub- 
grid models for SN and AGN feedback, and their influence on the 
galaxy population. Therefore, we will first describe the parameters 
and subgrid models used in the reference simulation. 

The simulation explicitly follows the 1 1 elements H, He, C, 
N, O, Ne, Mg, Si, S, Ca and Fe. Radiative cooling and heating 
are calculated element-by-elem ent in the presence of t he Cosmic 
Microwave Background and the lHaardt & Madaul d200ll) model for 
the UV/X-ra y background radiation from quasars and galaxies, as 
described in Iwiersma et ail d2009tj) . In these calculations, the gas 
is assumed to be in photo-ionization equilibrium and optically thin. 

In the centers of haloes the pressure is so high that the gas 
is expected to be in multiple phases, with cold and dense molec- 
ular clouds embedded in a warmer, more tenuous gas. This multi- 
phase structure is not resolved by our simulations (and the simu- 
lations lack the physics to describe these phases), so we impose 
a polytropic effective equation of state for particles with densi- 
ties « H > 10~' cirr 3 . These particles are also assumed to be star 
forming, as th is is the density required to form a molecular phase 
<Schavdl2004h . We set the pressure of these particles to P oc p^rff, 
where y cff is the polytropic index and p is the physical proper mass 
density of the gas. In order to prevent spurious fragmentation due 
to a lack of numerical resolution we set y eff = 4/3, as then the ratio 
of the Jeans length t o the SPH kernel and the Jean s mass are in- 
dependent of density dSchave & Dalla Vecchial2 008). The normal- 
ization of the polytropic equation of state is such that for atomic 
gas with primordial composition, the energy per unit mass corre- 
sponds to 10 4 K, namely (P/k = 1.08 x 10 3 K cirr 3 for n H = 10" 1 
cm" 3 ). Th e implementation of star forma tion is stochastic, as de- 
scribed in Sch ave & Dalla Vecchid (2008), with a pressure depen- 
dent SFR, obtained from local hydrost atic equilibrium and the ob- 
served Kennicutt-Schmidt law (KS-law Kennicutt 1998). 

For the same 1 1 elements that we use for the cooling, we fol- 
low the production by AGB stars and by Type la and Type II (in - 
cluding Type Ib,c) SNe, as described in Iwiersma et al. (2009b). 
The star parti cles are a s sumed to be simple stellar populations 
(SSPs) with a IChabrierl d2003l) initial mass function (IMF). SN 
feedback is implemented kinetically. After a short delay of 30 Myr, 
corresponding to the maximum lifetime of stars that end their lives 
as core-collapse SNe, newly formed star particles inject kinetic 
energy into their surroundings by kicking a fraction of their SPH 
neighbours in random directions. Each SPH neighbour i of a newly 
formed star particle j has a probability of J]nij/ ' m, of receiv- 
ing a kick with a velocity v w . Our reference model uses rj = 2 and 
v w = 600km/s, which for our assumed IMF corresponds to ~ 40% 
of the available energy from SNe being injected as winds. See Ta- 
bleQjfor the parameters that are varied in this paper. 



The only signicant discrepancy is in erg, which is 8 per cent, or 2.3cr, 
lower than the value favoured by the WMAP 7-year data 



2.2 Halo identification 

Haloes are identified using a Friends-of-Friends (FoF) algorithm, 
which links together all dark matter particles which are closer to 
each other than the linking parameter (b = 0.2 times the mean inter- 
particle distance). FoF ide ntifies iso-overdensity contours of S = 
(p-p)Jp = 3/(2nb 3 ) = 60 dDavis et alJl985l ; lLacev & Cofell994l) . 
Baryonic particles are linked to their nearest dark matter particle 
and belong to the same group, if any. 

Following the convergence tests presented in Appendix [A] we 
only include haloes that contain at least 100 star particles when 
considering halo properties as a function of stellar mass and we 
use a minimum of 2000 dark matter particles when we plot prop- 
erties against halo mass. These two cuts produce nearly identical 
halo samples in the reference simulation and ensure that only well 
resolved haloes are considered. In Appendix|B]we confirm that our 
results are robust to changes in the definition of halo mass used. 

Whenever we show the correlation between two halo prop- 
erties, the plot consists of lines that connect the medians of bins 
evenly spaced in the quantity plotted along the horizontal axis. Each 
bin contains at least 30 data points. If there are fewer than 30 points 
in a given bin, it is extended until it includes 30 objects. The last 
bin may contain between and 30 objects. We bin the data starting 
from the high-mass end. There, the difference in mass for two con- 
secutive haloes is much larger than at the low-mass end, and in this 
way we are sure that the value of the mass at the high-mass end of 
the plots is always the mean of the mass of the 15th and 16th most 
massive systems. 



2.3 Comparing simulations to observations 

In this section we describe how we compare our simulated galaxy 
population to observed stellar masses and SFRs. To convert ob- 
servationally inferred stellar masses and SFRs from the cos- 
mology assumed in the literature to our cosmology, we mul- 
tiply them by the square of the ratio of luminosity distances, 
[rft,ourcosm(z)/^L,obsco8m(z)] 2 . The subscripts 'our cosm' and 'obs 
cosm' denote our cosmology and the cosmology under which the 
observations are transformed into masses/SFRs, respectively. 

Note that we are using FoF halos, so all satellites are added to 
the central galaxy. We show in Appendix |B]that it makes very little 
difference to treat the satellites separately. 

We compare our simulated SFRs to those observed by 
iDaddi et al. (2007), who measured obscured and unobscured star 
formation by taking SFRs from both the UV and I R for /("-selected 
sBzK galaxies (star forming, see Da ddi et alj2004l) in the GOODS 
fields at z ~ 2. The median of the observed SFR as a function of 
stellar mass is well fit by SFR = 2 50 x (M./10 11 M B f s . The IMF 
assumed in the observations is the Salpeter (1955) IMF, whereas 
our stellar masses and SFRs are based on the Chabrier d2003l) IMF. 
We therefore divide the observationally inferred SFRs by a factor 
1.65, which is the asymptotic (reached after only 10 8 yr) ratio of 
the number of ioni zing photons per unit stellar mass predicted by 
iBruzual & Charlotl d2003l) for a constant SFR. 

For stellar masses, the IMF conversion factor is more sensi- 
tive to the age of the population and the observed rest-frame wave- 
length. As the light in most wavelength bands is dominated by mas- 
sive stars and the high-mass ends of both the Salpeter and Chabrier 
IMFs are power laws with very similar power law indices, we use 
the same factor of 1.65 as we used for the SFRs. For very old pop- 
ulations observed in red wavelength bands (tracing stellar continua 
rather than dust emission) the conversion factor should be differ- 
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Table 1. Overview of the simulations and the input physics that is varied in this paper. Bold face indicates departures from the reference model. The first 
column gives the name of the simulation; the second column denotes the type of SN feedback, either kinetic, thermal or none. The third column shows the 
wind velocity for kinetic feedback models (the circular velocity is defined as v c = ^GM yn /R y i r and the velocity dispersion <x is related to the gravitational 
potential (D: <x = V-<I>/2); the fourth column indicates, for kinetic feedback models, the wind mass loading. The fifth column indicates whether the winds 
are decoupled from the hydrodynamics, the sixth column indicates whether metal-dependent cooling is followed in the simulation and the seventh column 
indicates the stellar IMF(s). The eighth column shows which simulations include AGN feedback and the last column specifies the section in which each 
simulation is discussed. 



Name 


Kinetic/thermal 


^'wind 




Winds 


Z 


IMF 


AGN 


Sect. 




SN feedback 


(km/s) 


decoupled? 


Cooling 




feedback? 




REF 


Kinetic 


600 


2 


no 


yes 


Chabrier 


no 


All 


NOSN JVOZCOOL 


None 


n.a. 


n.a. 


n.a. 


no 


Chabrier 


no 




NOZCOOL 


Kinetic 


600 


2 


no 


no 


Chabrier 


no 


EH 


WML4 


Kinetic 


600 


4 


no 


yes 


Chabrier 


no 


El 


WML8V300 


Kinetic 


300 


8 


no 


yes 


Chabrier 


no 


EH 


WML4V424 


Kinetic 


424 


4 


no 


yes 


Chabrier 


no 


eh 


WML1V848 


Kinetic 


848 


1 


no 


yes 


Chabrier 


no 


EH 


WDENS 


Kinetic 


~ C s <x p 1 ' 6 


~ p-!/3 


no 


yes 


Chabrier 


no 


EH 


DBLIMFML14 


Kinetic 


600 


2, 14 


no 


yes 


Chabrier, top-heavy 


no 


EH 


DBLIMFV1618 


Kinetic 


600, 1618 


2 


no 


yes 


Chabrier, top-heavy 


no 


EH 


DBLIMFCONTSFV1 61 8 


Kinetic 


600, 1618 


2 


no 


yes 


Chabrier, top-heavy 


no 


EH 


WVCIRC 


Kinetic 


5v c /V2 


150/(V2v c ) 


no 


yes 


Chabrier 


no 


EH 


WPOTNOKICK 


Kinetic 


3cr 


150/(V2v c ) 


no 


yes 


Chabrier 


no 


EH 


WPOT 


Kinetic 


5<r 


150/(V2v c ) 


no 


yes 


Chabrier 


no 


EH 


WHYDRODEC 


Kinetic 


600 


2 


yes 


yes 


Chabrier 


no 


EH 


WTHERMAL 


Thermal 


n.a. 


n.a. 


n.a. 


yes 


Chabrier 


no 




AGN 


Kinetic 


600 


2 


no 


yes 


Chabrier 


yes 


EH 



ent. We verified that the K-band mass-to-light ratio is about a factor 
1.65 smaller for a Chabrier than for a Salpeter IMF for SSPs and 
constantly star forming pop ulations, for the full range of ages and 
metallicities available in the lBruzual & Charlotl 12003) population 
synthesis package. We therefore also divide by a factor of 1.65 to 
convert stellar masses from the Salpeter to the Chabrier IMF. 

We compa r e to the galaxy stellar mass function of 
iMarchesini et alj J2009h . a combined sample, using the deep near- 
infrared Multi-wavelength Survey by Yale-Chile, the Faint Infrared 
Extragalactic Survey and the Great Observatories Origins Deep 
Survey-Chandra Deep Field South su rveys at z ~ 2. Specific ally, 
we compare to the 1/V max results of Marchesini et alj ( l2009h . in- 
cluding all of their uncertainties, with the exception of bottom- 
light IMFs. The reason for this choice is that models with bottom- 
light IMFs dominate the systematic errors and represent a more 
extreme assumption than the variations in the other quantities. Ad- 
ditionally, theoretical justifications for using bottom-light IMFs ex- 
ist, thus far, only at high redshift (Dave 2008; van Dokkum 2008; 
IWilkins et all 12008). We obtain the error bars on the observed 
data points by considering all of the sources of random errors in- 
cluded in the observational data points, which include: Poisson er- 
rors on the number counts, cosmic variance and the random er- 
rors from the use of photometric redshifts. We add the sources 
of random error in quadrature and, in addition, linearly add the 
maximum of the system atic errors in the same mass bins, just as 
IMarchesini et alj J2009h . The systematic errors include the system- 



packages i 


Marchesini et al.l 


iMarastonl 


2005, Chariot & 



Bruzual, in prep.), varying the metallicities of the stellar popula- 
tions, jmdthe use of di fferent extinction curves (Milky Way f rom 
Allenlll976L SMC f rom lPrevot et alJll984llBouchet etai]|l985l and 



Calzett i et alfeOQOl) . 



Because our redshift of interest (z = 2) c orresponds to the 
bound ary between two of the redshift bins of Marchesini et al. 
(2009), which are 1.3 < z < 2 and 2 < z < 3, we weigh the 
averaging to the sizes of the redshift intervals (weights 1 .2 Gyr and 
0.8 Gyr respectively), which gives results consist ent with the z = 2 
results of the Newfirm Medium-Band Survey ( Mar chesini et al.1 
120101) . Additionally, the observed mass bins are not exactly the 
same size in both redshift intervals, although the differences are 
very small. Observed mass bins are of size 0.3 dex (at 1.3 < z < 2) 
and 0.29 dex (at 2 < z < 3). We interpolate the observed mass bins 
to a consta nt size of 0.3 dex . The r esulting 1/V max estimate of the 
analysis of Marchesini et al. (2009) is shown as the yellow shaded 
regions in the bottom-right panels of Fig.[T] and Figs.[3]-[9] labelled 
(I). 

March esini et~afl j2009l) use a diet Kroupa IMF, for which the 
correction factor to our Chabrier IMF is very small (~ 1.03, diet 
Kroupa being slightly more massive for the same flux). As this 
number is also derived from population synthesis packages, which 
come along with their own uncertainties, we chose not to convert 
masses for the small difference in IMFs. We do correct the masses 
for the difference in luminosity distances as described earlier. Num- 
ber densities also need to be converted, as the volume at a given 
redshift is different for different angular diameter and co-moving 
distances. Therefore, the number density (0„) is corrected for the 
ratio of volume elements (which is a function of the assumed cos- 
mology). 

When we compare our simulated galaxies to observations of 
the molecular gas mass in galaxie s, we use a sub-set of the compi- 
lation used in lGenzel et al. d2010h . As their total range of redshifts 
is very large, we chose to use only the two sub-sets directly above 
and be low z = 2. The near-IR long-slit Ho- sample from lErb et al.1 
(2006) contains 1 1 galaxies with a mean redsh ift of 2.3 and i s orig- 
inally drawn from the BX selected sample of ISteidefe t al. (2004) 
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and lReddvet~ai1d2005t) . The 4 galaxies with a mean re dshift of 1.5 
that we inclu de are star-forming BzK galaxies from iDaddi et al.l 
d2008l l201Ch . For both samples the molecular gas masses come 
from CO measurements with the Plateau de Bure interferometer 
and the stellar masses come from optical/UV SED fitting, under 
the assumption of a Chabrier IMF. We therefore only have to cor- 
rect for the difference in luminosity distance, as explained above. 
Our simulations do not explicitly calculate the molecular fraction 
in the high-density gas. Nevertheless, the galaxies we compare to 
are at the very highest stellar mass end of our simulated galaxies. 
For those galaxies the molecular gas fraction (by mass) is thought 
to be very high due to high pressures in the ISM, so we assume that 
all gas on the imposed equ ation of state is molecular. A recent study 
bv lNaravanan et aU d2012h suggests that with an improved conver- 
sion factor from CO to H2, gas f ractions of these gal axies are likely 
somewhat lower than quoted bv lGenzel et al.l d2010h . 



3 PHYSICAL PROPERTIES OF SIMULATED GALAXIES 
AS A FUNCTION OF DARK MATTER HALO MASS 

Fig. Q] shows, as a function of total halo mass the nine differ- 
ent galaxy properties we consider in this paper: medians of stel- 
lar mass fraction (/ star = M star /M tot , panel A), SFR (panel B), 
baryon fraction (fb myo „ = Mb Myo „/M tol , panel C), fraction of star- 
forming gas (/ism = MisM/Mtat, panel D) and gas mass fraction 



(/ga 



(M 



gas, total 



' A^isM)/Mot> panel E). Then, as a function of 



stellar mass: medians of the molecular gas mass in the ISM (panel 
F), specific SFR (sSFR = SFR/M,, panel G), the inverse of the gas 
consumption timescale (SFR/M ISM , panel H) and galaxy number 
density (the galaxy stellar mass function, panel I). 

In each panel of Fig. [T] we show the properties of the galax- 
ies in the reference simulation with a black curve. The blue curves 
show results from the simulations analyzed in Paper II, which in- 
clude changes in the cosmology, ISM physics, star-formation law 
and IMF. The red curves show the simulations discussed in this pa- 
per (see section [4j, which comprise changes in the feedback pro- 
cesses that are modeled. The red and blue lines are included in 
Fig. E to give a visual impression of the magnitude of the effects 
that we are considering. Here, we note that the magnitude of the 
differences between simulated galaxy properties due to uncertain- 
ties in the sub-grid physics are larger than the differences that arise 
from changin g the hydrodyn amics method from SPH to a (moving- 
) mesh code dSpr ingel 2010i), as e.g. shown by Vogelsberger et al.l 
d2012h : lKeres et al.l d2012Tu nd lScannapieco et alj d2012l) . The rela- 
tive accuracy of different numerical techniques can nevertheless be 
important, but an investigation of this issue falls outside the scope 
of this paper. 

In the remainder of this section we describe and explain the 
trends seen in the properties of the galaxies in the reference simu- 
lation, and comment briefly on the effects of the physics variations. 
Note that the first five panels (galaxy properties as a function of DM 
halo mass) and the last four panels (galaxy properties as a function 
of stellar mass) are resolved down to slightly different resolution 
limits, as discussed in Appendix lAl 



3.1 Properties as a function of halo mass 

In panel (A) of Fig.[T]we plot the stellar mass fraction as a function 
of total halo mass. Stellar mass and total mass are tightly (and al- 
most linearly) correlated so differences between the models are em- 
phasized by plotting the ratio of the masses. In the reference model, 



the fraction of the total mass locked into stars increases smoothly as 
a function of halo mass, from ~1% in the smallest resolved haloes 
(log 10 (M lo ,/M o ) ~ 10.5), to 4% in the highest mass haloes. It is 
difficult to make a quantitative comparison with gala xy formatio n 
efficiencies quoted in much of the literature (e.g. lGuo et alj|201oh . 
as the definition of halo mass used in this study is different from 
the spherical overdensity halo definitions used in these works, but 
the qualitative trend is the same. As we move towards higher halo 
masses, a fractionally larger proportion of the mass is locked into 
stars. It is clear from inspection of panel (A) that the changes in 
physics considered in this paper (red curves) have a large effect on 
the Mh a io - M, relation (much bigger than the simulations discussed 
in Paper II; blue curves). We discuss the reasons for this sensitivity 
in Sec. [4] The one simulation that is a strong outlier, with much 
higher stellar mass at fixed halo mass, is the simulation that ne- 
glects SN feedback and, as such, contains no mechanism to prevent 
the catastrophic cooling of gas into stars. 

In panel (B) we show the integrated SFR inside haloes as 
a function of their total mass. In all simulations the SFR is an 
increasing function of halo mass in the mass range probed here 
(10.5 < log 10 (M lo ,/M o ) < 12.0). As in panel (A), the vast major- 
ity of the physics variations considered in Paper II have very little 
effect, but the variations discussed here are very important. This im- 
plies that the sub-grid physics employed to model the ISM, other 
than that related to feedback, does not affect the large-scale, stel- 
lar properties of the galaxies. The outlier among the red curves in 
this panel is again the simulation that neglects SN feedback, which 
shows that feedback is an important component of the sub-grid 
models. 

Panel (C) shows the baryon fractions of the haloes as a func- 
tion of halo mass. The universal baryon fraction appropriate for 
our default cosmology (O b /f2 m = 0.18) is over-plotted as a hori- 
zontal, dashed line. In the default simulation, feedback processes 
are able to eject baryons from the halo very easily in low-mass 
haloes, leading to low baryon fractions (0.6 dex below universal) 
for M (ot ~ 10 10 5 M o , whereas they become gradually less efficient 
as the halo mass increases. Again, the physics changes presented 
in this paper have a much more significant effect on this relation 
than those discussed in Paper II. Panels (D) and (E) show the frac- 
tion of the mass that is in gas in the ISM (which is sensitive to the 
physics variations in Paper II) and the rest of the gas in the haloes, 
respectively. In general, both are increasing functions of the total 
mass. 



3.2 Properties as a function of stellar mass 

The mass in the ISM as a function of the galaxy stellar mass is 
slightly sub-linear, as depicted in panel (F). This reflects the fact 
that at higher masses, gas is more efficiently transformed into stars. 
Our simulations slightly under-predict the molecular gas mass as a 
function of stellar mass when compare d to the sub-sample of the 
galaxies studied bv lGenzel etai]d20ld) . 

Both the stellar mass (panel A) and SFR (panel B) of galax- 
ies are correlated almost linearly with halo mass, suggesting, that 
there is a linear relation between the stellar mass and the SFR. This 
linear trend masks important differences between the simulations, 
so to examine galaxy properties in more detail, we use the galaxy 
specific SFR (sSFR=SFR/M»). 

Panel (G) shows the median sSFR as a function of galaxy stel- 
lar mas s. In this panel, the black, dotted line shows the observa- 
tions of IDaddi et ail (2007). The scatter in the data is not shown, 
but is constant at approximately 0.2 dex. The scatter in the simu- 
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Figure 1. Median relations between halo properties in all the simulations described in this work (red lines) and in Paper II (blue lines). The reference model 
is shown as a black curve in each panel. In Sec. 4 we consider subsets of simulations in more detail. In the first five panels we show, as a function of total halo 
mass, tmedians of stellar mass fraction (/ star = M st ar/M 10t , panel A), SFR (panel B), baryon fraction (/baryon = ^bai-yon/Mot. panel C), fraction of star-forming 
gas (/ism = ^ISM/Mot, panel D) and gas mass fraction (/ gas , halo = ('Wgas, total _ ^lSM)/A/tot: panel E). The last four panels show, as a function of stellar mass, 
medians of the molecular gas mass in the ISM (panel F), specific SFR (sSFR = SFR/M,, panel G), the inverse of the gas consumption timescale (SFR/M[sm, 
panel H) and galaxy number density (the galaxy stellar mass function, panel I). As described in the text, we show medians in bins along the horizontal axes for 
all haloes that satisfy the convergence criteria that apply to that specific panel. The horizontal, dashed line in panel (C) shows the universal baryon fraction for 
our chosen cosmology, the data points in panel (F) show a sub-set of the compilation studied in Genzel et al. (2010), the dotted black line in panel (G) shows 
the stellar mass - sSFR relation from the GOODS field (Daddi et al. 2007) and the shaded yellow region in panel (I) shows the galaxy stellar mass function of 
Marchesini et al. (2009). 
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lation data points is similar, although somewhat smaller (~ 0.1 to 
0.15 dex). We plot the observed relation over the mass range that is 
observed, and do not extrapolate beyond there. 

The trend seen in the simulations is such that at low masses the 
sSFR is an increasing function of mass, but at ~ 10 10 M Q the trend 
reverses and the sSFR decreases with increasing stellar mass. At the 
highest masses the relation is declining more steeply than observed. 
As discussed in Sec. 13. II most of the physics variations discussed 
in this paper have little effect on the integrated stellar properties of 
the haloes. The slight underestimate of the SFR is likely connected 
to the underestimate of molecular mass shown in panel (F). 

Because SFRs are expected to be more strongly influenced by 
the mass of available star-forming gas than by the amount of stars 
already formed, we define a second normalized SFR, the ratio of 
the galaxy SFR to its star-forming gas mass. This is the inverse of 
the time needed to convert the present reservoir of star forming gas 
into stars, assuming that the present SFR is maintained. We plot 
the inverse of the gas consumption time scale as a function of stel- 
lar mass in panel (H). The gas consumption time is a monotonically 
decreasing function of stellar mass in the simulations. In the cores 
of the more massive haloes, gas reaches higher densities and pres- 
sure, so gas is converted more efficiently into stars. 

In panel (I) we examine the galaxy stellar mass function: the 
number density of galaxies as a function of their stellar mass. Th e 
yellow, shaded regions show the data from lMarchesini et alj J2009I) . 
When the observational uncertainties are taken into account, most 
of the simulated mass functions fall well within the observed range. 
We note that, at low masses, our simulated stellar mass func- 
tion is steeper than most faint-end slopes of de rived Schechter 
funct i on parametrizations at simil ar redshifts (e.g. lConselice et all 
l2005MMarchesini et alj2009ll2010l) . but that this discrepancy exists 
largely outside the range of stellar masses where the mass function 
has been observed. Our simulation box size is too small to form the 
rare objects that populate the exponential cutoff of the stellar mass 
function. 

Comparing the behaviour of the reference simulation between 
panels (G) and (I) reveals an interesting behaviour: although the 
galaxy sSFR is a factor of a few too low at z = 2, the sim- 
ulations form enough galaxies of all masses, and possibly too 
many low-mass systems (a similar behaviour was also noted by 
IChoi & Nagaminell2012h . This may indicate that the observations 
are not internally consistent or that other physics variations are re- 
quired. 



4 ISOLATING THE EFFECTS OF THE INPUT PHYSICS 

In this section we discuss each of the variations of the input physics. 
Table Q] summarizes the simulations discussed in this paper and in- 
dicates in which subsection each one is discussed. Bold face values 
indicate differences between each simulation and 'REF' . In this pa- 
per we discuss how metal-line cooling and the inclusion of kinetic 
SN feedback affect galaxy properties (Sec. 14. it ; how the parame- 
ters of the kinetic SN feedback model affect the results fSec. l4.2l l: 
the effect of a top-heavy IMF at high pressures (Sec. 14. 3b ; winds 
with 'momentum-driven' scalings that depends on galaxy proper- 
ties (Sec. l4.4t : the effects of decoupling the winds from the hydro- 
dynamics (Sec. \4.5h the effect of injecting SN energy thermally 
(Sec.|4.6ll: and the effects of strong, AGN feedback (Sec. |4.7> . 

A graphical representation of the gas density of a galaxy 
formed in a representative set of models is shown in Fig. [2] The 



galaxy resides in a halo of total mass 



10'- 



Mra. It was first 



identified in the 'REF' simulation, where its position (defined as 
the centre of mass of all particles within 10% of the virial radius) 
is determined. The line of sight is along the z-axis, which is al- 
most perfectly aligned with the angular momentum vector of the 
gas within 10% of the virial radius (cos(0) = 0.994). For the other 
simulations the image is centered on the same position, showing 
the remarkable similarity in the positions and orientations of the 
galaxies. It is immediately clear that changing the feedback model 
can have a large effect on the galaxy morphology. Going so far as 
almost completely destroying the galaxy disk in the strongest cases 
('AGN' and 'DBLIMFV1618'). 

4.1 Metal-line cooling 

In this set of simulations we investigate how gas cooling through 
metal lines affects the galaxy population. Metal-line cooling is the 
dominant mechanism by which enriche d gas can cool in the te m- 
perature range 10 5 K < T < 10 6 K (e.g. IWiersma et alj|2009"3) . If 
gas shock heats to high temperatures while accreting onto galax- 
ies, neglecting metal-line cooling will greatly decrease the sup- 
ply of gas that can co ol out of haloes to fuel the galaxy (see 
van de VoortetaT]|201l[ for a comprehensive discussion of galaxy 
fueling in these simulations). In order to isolate these effects we 
compare the 'REF' simulation to one in which cooling through 
metal-lines was switched off ('NOZCOOL'), and to a further simu- 
lation in which both cooling through metal-lines and SN feedback 
were switched off ('NOZCOOLJVOSN'). The morphology of a typ- 
ical massive galaxy in these simulations can be seen in Fig.[2] Turn- 
ing off metal-line cooling reduces the extent of the gaseous disk in 
this massive system as less gas cools out of the halo into the galaxy. 
Neglecting SN feedback leads to a hugely centrally concentrated 
galaxy as there is no mechanism to eject low-angular momentum 
material or to stop the fragmentation, and associated angular mo- 
mentum losses. 

To isolate the effect of metal-line cooling, we can compare 
the 'REF' simulation (solid, black curve) to the 'NOZCOOL' sim- 
ulation (dotted, red curve) in Fig. [3] Panel (B) of Fig. [3] shows 
the effect of metal-line cooling on the SFRs of galaxies. In gen- 
eral, metal-line cooling increases SFRs, because cooling rates in- 
crease with increasing metallicities. The magnitude of the dif- 
ferences between the simulations increases with halo mass be- 
cause both the fracti on of gas that is shock h eated when it en- 
ters the galaxy (e. g. [Birnboim & Dekell 120031 : iKeres etalj 120051: 
lOcvirk et all lioO^Tvan de Voort et al.1 1201 ll) and the halo virial 
temperature increase with mass, making the effect of metal-cooling 
more pronounced. The same trends are evident in other gas proper- 
ties as a function of halo mass and stellar mass. Without metal-line 
cooling, less gas is able to cool onto the galaxy (panel D), and hence 
stellar masses (panel A) are lower. 

As a function of stellar mass, we see that sSFRs are higher 
(panel G) and gas consumption timescales are shorter (panel H) in 
the 'REF' than in the 'NOZCOOL simulation, as gas can cool more 
quickly into the galaxy when metal-line cooling is allowed, leading 
to higher gas densities and more efficient star formation. The stellar 
mass functions of the two simulations (panel I) are almost identical 
up to M star = lO 9 5 M , but above this mass the 'NOZCOOL' sim- 
ulation lies systematically 0.1-0.2 dex below 'REF'. This occurs 
because the normalization of the M stal - M tot relation is decreased 
in 'NOZCOOL' , and so a given stellar mass corresponds to a larger, 
rarer halo. 

Turning our attention now to the effect of SN feedback by 
comparing the 'NOZCOOL' simulation (dotted, red curve) to the 
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Figure 2. A graphical representation of a galaxy in a halo of 10 12 ' 5 M© in 15 of our simulations at redshift 2. The colour coding denotes the gas density in a 
slice of 100 ft -1 kpc thickness, divided by the mean density of the universe. All frames are 100 co-moving /i~'kpc on a side and are centered on the position 
of the galaxy in the 'REF' simulation. All frames have a thickness of 100 co-moving /i~'kpc. The orientation of the line of sight is along the z-axis, which is 
almost perfectly aligned with the angular momentum vector of all material inside 10% of the virial radius of this galaxy in the 'REF' simulation. It is clear 
that different feedback prescriptions can have a large effect on the size of the galactic disk. Note that 'WML4' is not included here. That simulation is also 
discussed in Paper II, and it is included there. 



'NOZCOOLJJOSN' simulation (dashed, blue curve), we see that 
neglecting SN feedback leads to dramatic changes in the galaxy 
population as with no energy injection there is no process that can 
prevent gas from cooling directly into the galaxy. 

Notably, panel (C) demonstrates that, unlike all of the other 
simulations we consider, 'NOZCOOLJVOSN' has baryon fractions 
above the universal value. This is due to the fact that star forma- 



tion and cooling allow more gas to be pulled into the potential 
well of the galaxy. Non-radiative simulati ons predict bary on frac- 
tions of ~ 0.8 times the universal value dCrain et alj|200l but see 
iKravtsov et al ]|2005l who find a value much closer to one in a non- 
radiative mesh-based simulation), as pressure support forces more 
gas outside of haloes. 

Interestingly, we note that (unlike all of the other simulations 
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Figure 3. As Fig.[T] but showing only the subset of simulations in which the metal-line cooling and/or kinetic SN feedback are turned off. The solid, black 
curve shows the 'REF' simulation. The red, dotted line shows the effect of turning off only metal-line cooling ('NOZCOOL'). The blue, dashed curve shows 
the effect of switching off both metal-line cooling and SN feedback ('NOZCOOLJJOSN'). The effect of the SN feedback can thus be isolated by comparing 
the blue, dashed and red, dotted curves. 



considered here), the slope of the M slal - sSFR relation (panel G) 
in the 'NOZCOOLJVOSN' has an almost identical slope to that ob- 
served, although with a normalization that is lower by a factor of 
approximately two. That the slope of the observed M stllr - sSFR re- 
lation is so similar to that of the simulation without feedback may 
indicate that perhaps feedback is ineffective at high galaxy masses, 
as is the case for the reference simulation above 10 10 M o (see also 
Paper II). The similarity to the observed relation is striking. The 
galaxies in the observations are selected to be star forming, so one 



would naively expect the star formation in those galaxies to be reg- 
ulated by feedback. 

4.2 Winds with constant energy per unit stellar mass formed 

Winds driven by feedback from star formation have become a fun- 
damental ingredient of the galaxy formation picture by preventing 
too much gas from being locked into stars. In the 'REF' simula- 
tion, such winds are assumed to be driven by SNe and are imple- 
mented kinetically. That is, the effect of SNe is to 'kick' some 
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Figure 4. As Fig. [T] but showing only a set of simulations in which the initial wind velocity and the initial wind mass loading are varied. The reference 
simulation (black, solid curve) has a mass loading of i] = 2 and a wind velocity of 600 km s . The simulations shown by the red-dotted, blue-dashed and 
green-dot-dashed lines show simulations that use the same total energy per unit stellar mass but different combinations of the mass loading and wind velocity. 
The simulations shown have mass-loadings of 1 ("WML1V848' ; red, dotted curve), 4 ('WML4V424' ; blue, dashed curve) and 8 CWML8V300' ; green, dot- 
dashed curve). The magenta dot-dot-dot-dashed line represents a simulation ('WDENS') which has a mass loading and velocity dependent on the local density, 
such that the energy in the wind is still the same and the initial velocity is proportional to the local sound speed. Changing the parameters of the kinetic 
feedback model can change the galaxy stellar properties by up to an order of magnitude, even when the energy that is injected is kept constant. The halo mass 
above which SN feedback becomes inefficient at removing gas from the galaxy (as seen by a steeper rise in the galaxy SFR with halo mass) depends primarily 
on the SN-driven wind velocity. 



amount of gas. This model contains two free parameters: the ini- 
tial mass-loading 77 = M win( |/M, and the initial wind velocity v w , 
which remain essentially unconstrained. We therefore compare a 
series of four simulations that use the same SN energy per unit 
stellar mass formed (i.e. rjv^ is the same as for 'REF'), but dis- 



tribute it differently between the velocity and mass loading. The 
mass loadings in the four simulations are 1, 2 (i.e. l REF'), 4 and 8, 
with corresponding velocities of 848, 600, 424 and 300 km/s, re- 
spectively. The parameters of the wind model are contained in the 
simulation name. For example, the simulation 'WML1V848' uses 
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a mass-loading of 1 and a wind velocity of 848 km s -1 . A fifth 
simulation with kinetic feedback with the same feedback energy, 
' WDENS' , uses a wind velocity that scales in proportion to the lo- 
cal gas sound-speed (oc p 1/6 ) and a mass loading that scales such 
that the total injected energy is constant (oc p~ 1/3 ). The normaliza- 
tion of the parameters in the ' WDENS' simulation is such that the 
wind velocity and mass loading are the same as in the reference 
model if the gas density equals the star formation threshold, i.e. 
«h =0.1 cm 4 . For illustrative purposes, we also consider a simu- 
lation in which the available SN energy relative to 'REE' is doubled 
(by doubling the mass-loading to 4 and keeping the wind velocity 
the same). This simulation is termed 'WML4' . 

We first compare the set of simulations in which we hold the 
amount of SN energy that is injected per unit stellar mass constant, 
but vary how the energy is distributed between the initial mass load- 
ing and wind velocity. We show in Fig.|4]that the simulation with 
the highest wind velocity ( ' WML1 V848; red, dotted curve) sup- 
presses SF much more effectively in high-mass objects than the 
simulations with lower wind velocities (e.g. 'WME8V300'). How- 
ever, at low galaxy masses, this trend reverses and it is the simu- 
lations with the lowest wind velocities that suppress SF most ef- 
ficiently (panels A, B and D). At the low-mass end, the amount 
of suppression is set by the amount of gas that is kicked, where 
models with large mass-loadings are able to inject large amounts 
of gas into the wind, much of which escapes into the halo or be- 
yond, leading to these models being most effective. However, in 
massive haloes a large wind velocity is required to drive gas out of 
the galaxy, so the models with large mass-loadings become inca- 
pable of driving winds and form more stars in high-mass objects. 
There exists a clear transition mass above which, for a given wind 
velocity, SN feedback is no longer able to drive gas from the galaxy 
(see also lSpringel & Hernqu ist 2003b). This is visible both in the 
properties of the galaxies (panels F-H), and the properties of the 
galaxy haloes (panels A-E), where above a mass (that depends on 
the initial wind velocity) the results tend towards those in a simula- 
tion that includes no feedback (Fig.[3j blue, dashed curve). 

iDalla Vecchia & Schavd |2008) showed that the reason why 
effective feedback in more massive galaxies requires higher veloc- 
ity winds is not directly related to the requirement that the winds 
overcome the gravitational potential. Instead, they found that the 
outflows are already quenched in the ISM and only if gas pressure 
forces are taken into account. Because the pressure in the ISM in- 
creases with the depth of the potential well, more massive galaxies 
require hig her velocity winds. A deeper understanding was pro- 
vided by IDalla Vecchia & Schavel d2012h . who demonstrated ana- 
lytically, and confirmed numerically, that the outflows are quenched 
due to radiative losses in the shock-heated ISM. Higher velocities 
imply higher post-shock temperatures and hence longer cooling 
times. As the ISM of more massive galaxies is denser, higher veloc- 
ities are required for effective feedback in more massive galaxies. 
They pointed out that cosmological simulations overestimate the 
radiative losses due to their finite resolution and showed that mod- 
els with different wind parameters, but the same amount of energy 
per unit stellar mass, make converging predictions if the resolution 
is sufficiently high. The resolution of our simulations is low com- 
pa red to the analytic estimates of the required resolution provided 
bv lDalla Vecchia & Schavel d2012h . which explains our finding that 
galaxy properties are sensitive to the choice of initial wind velocity. 

A second instructive comparison is between 'REE' and the 
simulation that injects twice as much SN energy per unit stellar 
mass, 'WML4' (orange, long-dashed curves). This simulation uses 
the same wind velocity as 'REE' and so the SN-driven winds be- 



come ineffective at the same mass. This can be seen by comparing 
'REE' (black, solid curve) to WML4 in panel (B), where above 
M tot = 10"' 25 M the two simulations converge towards each other 
with increasing mass, but below this mass, 'WML4' can suppress 
SF significantly more efficiently than 'REF' . We can see this more 
quantitatively in panel (A), where we see that at low masses (where 
the feedback is effective), / slar is almost exactly a factor of two 
lower than in the 'REE' simulation. This implies that the star for- 
mation is self-regulating: a galaxy's SFR increases until the rate of 
energy injection into the ISM reaches a critical value. This critical 
rate of energy injection depends on the halo mas s and is presumably 
the rate for which outflows balan ce inflows (e.g. Schave et al ■l20ld 
iDave et ai] |2012l ; see also lBooth & Schavejl20ld for an analogous 
discussion of self-regulated black hole growth). Because the out- 
flow rate depends not only on t he rate of energy injection, but also 
on the initial wind velocity (e.g. lDalla Vecchi a & Schave 2008|), the 
same will be true for the critical SFR. Doubling the SN energy that 
is injected per unit stellar mass, while keeping the initial wind ve- 
locity constant, will half the SFR required to produce the same out- 
flow rate and will thus cause the galaxy to form half as many stars. 
This represents one of the fundamental conclusions of this paper, 
star formation is regulated by the interplay between the available 
fuel supply and feedback processes. 

Because a given galaxy requires a fixed rate of energy injec- 
tion to achieve a balance between inflows and outflows, increasing 
the efficiency of SN feedback causes the galaxies to decrease their 
gas reservoirs (panel D). Because the star formation law is non- 
linear and because the ISM of more massive galaxies tends to be 
denser, the decrease in the ISM is larger for lower mass galaxies. 

The gas consumption timescales of all models are very similar 
in the regime where SN feedback is effective (panel H), and at the 
point where feedback becomes incapable of lifting gas out of the 
galaxy, the gas in the galaxy is used up very efficiently. 

For the 'WDENS' model, the energy injected in the wind per 
unit stellar mass is also the same as in the reference model, but the 
initial wind velocity scales with the local sound speed. The relation 
between halo mass and SFR is even shallower than it is for the 
run with v„. = 848 km s , indicating that the feedback is effective 
for all haloes. At high masses, this is our most effective wind model 
with constant energy. This indicates that choosing the mass-loading 
and wind velocity carefully can make SN-driven winds capable of 
driving winds in all objects resolved in these simulations. 

As is clear from panel I of Fig. [4] decreasing the slope of 
the low-mass end of the stellar mass function can be attained by 
increasing the mass loading factor in constant energy winds. The 
highest mass loading still gives a low-mass end slope that is some- 
what steeper than power law fits to the low-mass end in the ob- 
servations, although the discrepancy only occurs on masses lower 
than those observed. Such efficient feedback in low-mass objects 
also boosts the sSFRs in the higher-mass objects for which feed- 
back from star formation is just becoming inefficient, as appears to 
be required by the observations (panel G). 

These simulations suggest that the fit to the observed mass 
function could be much improved by varying the mass load- 
ing factor with galax y mass (see also lOppenheimer et alj [20101: 
Okam oto et a ] l20ld : iBower etalJ |2012| ; IPuchwein & Springe! 



2012b . 



4.3 Simulations with a top-heavy IMF at high pressures 

Determining the stellar IMF represents a significant challenge out- 
side of the solar neighbourhood. The IMF determined locally is 
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Figure 5. As Fig. [T] but showing only the set of simulations in which a top-heavy IMF is used for star formation at high pressure. The 'REF' simulation 
(without top-heavy IMF) is shown by the black, solid curve. The extra SN energy available from a top-heavy IMF can either go towards increasing the mass 
loading ('DBLIMFML14' ; red, dotted curve) or the initial velocity of the winds ('DBLIMFV1618' ; blue, dashed curve). In the presence of a sudden change in 
the IMF at some pressure, one can either allow the star formation rate or the rate of formation of massive stars (which is what is observed) to be a continuous 
function of the pressure ('DBLIMFCONSTSFV1618' ; green, dot-dashed curve). Putting the extra SN energy into the mass-loading ('DBLIMFML14') does not 
increase the level of suppression of the SFR, whereas increasing the wind velocity does allow more gas to escape from galaxies. 



usually taken to be universal and is applied to all galaxies, at all 
redshifts. There is, however, some evidence suggesting that the 
IMF may be 'top-heavy' in extreme regions such as starbursts (e.g . 
iBaugh et al.l 120051 : lHabergham etal] l2010fc IWeidner et all 1201 lb . 
We performed a series of runs in which a different IMF is assumed 
in stars that formed in high-pressure environments. 

In these simulations, stars form with an IMF N/M oc M _1 
(compared with M~ 23 for the high-mass end of the Chabrier IMF) 



if the gas pressure exceeds P/k = 2.0 x 10 6 cur 3 K (evaluated at 
the resolution limit of the simulations). This process mimics the 
observation that stars that form in extreme environments (such as 
starbursts) appear to have an IM F that is flatter than t he Chabrier 
IMF we assume in 'REF ' (e.g. iMcCradv et alj[2003l ; IStolte et ail 
l2005l : lManess et alj|2007b . 

For a top-heavy IMF, the available energy from SNe per unit 
stellar mass formed is higher by a factor of ~ 7 than for the Chabrier 
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IMF. Numerically, we have some freedom in how this energy is 
distributed. It can be used either to increase the wind mass loading 
or the wind velocity. We therefore ran two additional simulations 
with a top-heavy IMF at high pressures. In the first simulation the 
wind velocity remains fixed at 600 km s _1 (the 'REF' value), but 
the mass loading is increased to ^ = 14 {'DBL1MFML14') for stars 
that form in high-pressure environments. In a second simulation, 
the mass loading is kept fixed at 77 = 2 (the 'REF' value) and the 
wind velocity is increased to v„, = 1618 km s -1 ('DBLIMFV1618') 
for stars that form at high pressures. 

When changing the IMF suddenly as a function of the pres- 
sure, it is not immediately clear how to treat the star formation 
law. The Kennicutt-Schmidt law is inferred from observations that 
probe only massive stars. The actual SFR, therefore, depends on 
how many low-mass stars are formed per unit mass of massive 
stars. When changing the IMF, the star formation law can be 
changed in two ways: 

(i) From observations there is no indication of a discontinuity 
in the formation rate of massive stars with pressure. Although this 
is most likely the result of the IMF being a continuous function of 
SFR or pressure (if there is a relation at all), we nevertheless im- 
plemented a model that changes the normalization of the KS-law 
such that the formation of massive stars is continuous, resulting in 
a discontinuous SFR as a function of pressure (the KS-law normal- 
ization drops at the pressure above which the IMF is top-heavy). 
This is the procedure we follow in most of our simulations, and is 
the assumption made by both the simulations 'DBL1MFML14' and 
'DBLIMFV1618'. 

(ii) If the (total) SFR as a function of pressure is continuous, 
the formation of massive stars must be discontinuous, given that 
we assumed the IMF to change suddenly above some critical pres- 
sure. We run one additional model under this assumption, termed 
'DBLIMFCONTSFV1618', which assumes that the total SFR as a 
function of pressure is continuous and the extra energy injected by 
a top-heavy IMF goes into increasing the wind velocity relative to 
the 'REF' simulation. 

When comparing simulations to observations, we do not cor- 
rect the stellar mass of simulations with a double IMF. On average, 
only ~ 10% of the star particles in the simulation box formed with 
a top-heavy IMF (this depends slightly on resolution and weakly 
on whether the rate of for mation of mas s ive st ars is a continuous 
function of the pressure). In Sch ave et alJd20l6T) we showed that at 
late times, this correction can be significant, but that at z = 2 the 
integrated SFR of the universe is the same, regardless of whether 
or not the SFRs of particles at pressures higher than the threshold 
pressure for the top-heavy IMF are corrected for another assumed 
IMF. Also, we demonstrate in Paper II, by considering simulations 
in which the normalization and slope of the star-formation law are 
varied, that the form of the assumed SF law does not strongly influ- 
ence the galaxy mass function or the SFR distributions in galaxies 
(it mainly affects the ISM fraction). At z = 2 we therefore expect 
any differences between the runs that use a top-heavy IMF and the 
'REF' simulation to be mostly due to the extra energy input from 
SN feedback and/or the increased rate of production of metals that 
results from a top-heavy IMF. 

We show results from these simulations in Fig. [5] In all 
of the plots there are two 'families' of curves. The simulations 
'REF' (black, solid curves) and 'DBLIMFML14' are, in almost all 
plots, similar to one another, and the simulations 'DBLIMFV1618' 
(blue, dashed curves) and 'DBLIMVCONTSFV1618' (green, dot- 
dashed curves) are also very similar to one another. The fact that 



'DBEIMFV1618' and 'DBEIMVCONTSFV1 618' are always very 
similar tells us that, if feedback is effective, whether or not we have 
a continuous SFR at the threshold pressure is only of minor impor- 
tance when it comes to galaxy properties, and as such we will not 
consider 'DBE1MVCONTSFV1618' any further. 

Next we compare 'REF' and 'DBLIMFML14' . As discussed 
in Sec. 14.21 SNe are only capable of driving winds out of a galaxy 
when they are launched at sufficiently high velocities. Therefore, 
'DBLIMFML14' behaves similarly to 'REF' because the two mod- 
els are capable of driving winds out of the same objects. It is appar- 
ent that 'DBLIMFML14' does form fewer stars than 'REF' at high 
halo masses (M tot > 10 M G ; panels A and B), even though the 
halo baryon fractions of the two simulations are very similar (pan- 
els C and E). That the two models differ only at high halo masses 
can be explained by the fact that the pressure in the ISM increases 
with the depth of the potential. Hence, only in massive haloes does 
a significant fraction of the ISM have pressures above the critical 
value above which we assume a top-heavy IMF. 

The reason why 'DBL1MFML14' has lower stellar mass frac- 
tion in massive haloes than 'REF' can be seen by comparing the 
gas consumption timescale (panel H) and the mass of gas in the 
ISM (panels D and F) for the same two simulations. Although there 
is significantly more gas in the ISM in 'DBLIMFML14' than in 
'REF', the gas consumption timescale is much longer. This may 
reflect the drop in the normalisation of the KS-law in high-pressure 
gas in 'DBLIMFML14' . Because the gas consumption time scale 
becomes extremely long in the high-pressure gas (~ 10 10 yr) and 
because the initial wind velocity is too low to be effective, self- 
regulation may be difficult to achieve. 

Finally, we compare 'REF' and the top-heavy IMF simulation 
in which the initial wind velocity is increased ('DBLIMFV1618'). 
It is immediately obvious that the effect of the top-heavy IMF is 
to strongly suppress star formation (panels A, B, F, G and H) as 
well as to effectively remove gas from galaxies (panel D) and from 
their host haloes (panels C and E). The effects increase with the 
mass, because the ISM pressure, and thus the fraction of stars form- 
ing with a top-heavy IMF, increase with the depth of the potential 
well. The effect of the top-heavy IMF is very similar to increas- 
ing the wind velocity (see 34. 2\ . In both simulations the SN-driven 
winds are able to drive out gas from all but the most massive ob- 
jects, but the stellar fractions are lower for 'DBLIMFV1618' than 
for 'WML1V848' . This is a result of the increase in the available 
SN energy per unit stellar mass, so a smaller star formation rate 
is required to drive winds from a galaxy and to effectively regu- 
late star formation. Finally, we note that the simulations that in- 
clude strong feedback from a top-heavy IMF suppress the ampli- 
tude of the galaxy stellar mass function below the observed points 
(panel (I)). This behaviour is seen in the other models that contain 
very strong feedback (e.g. 'WPOT' and 'WVCIRC in Sec.|4~4land 
'AGN' in Sec.[4?7t. 

We conclude that the effect of a top-heavy IMF at high pres- 
sures mainly reflects the increased efficiency of winds driven by 
massive stars. The extra energy available for driving winds is more 
important than the increase in the metal yields associated with a 
top-heavy IMF. 

4.4 'Momentum-driven' wind models 

Galactic winds can be driven by radiatio n pressure on dust gr ains in 
the wind, which drag along the gas (e.g. lMurrav e t al. 2005). Here 
the driving force of the wind is the radiation pressure which injects 
momentum into the outflow. If the cooling time in the shock-heated 
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Figure 6. As Fig. [T] but showing only the reference simulation ('REF'; black, solid curve) and the set of simulations in which 'momentum-driven wind' 
scalings are employed. In 'WVCIRC (red, dotted curve), the initial wind velocity scales with on the circular velocity of the halo the wind is launched from, 
while in 'WPOTNOKICK' and 'WPOT' the wind velocity is set by the the local gravitational potential (without and with an extra kick respectively, shown 
by the blue, dashed and green, dot-dashed curves respectively). In these models the energy injected into the wind per unit stellar mass is not constant and 
generally exceeds the energy in the reference simulation, as is evident from the increased suppression of star formation. 



gas is sufficiently short, then the wind will conserve its momentum 
and the wind is said to be 'momentum-driven'. Although cosmo- 
logical, hydrodynamical simulations have not yet tried to model 
this process directly by including radiation transport, simulations 
have been run that include the same type of kinetic feedback as 
used for SN-driven winds, but in which the initial wind velocity 
varies with local physical conditions or galaxy properties, and the 
initial mass loading varies as rj oc v w ' , such that the rate with which 
momentum is injected per unit stellar mass, oc r]v w , is constant. 



This scaling should be contrasted with that used for the 'constant 
energy' winds considered earlier, for which rj oc y^ 2 . 

We implemented several 'mo mentum-driven wind' models 
that are similar to those used by lOppenheimer & Pavel J2006| . 
2008). Here, the wind parameters var y with either the local poten- 
tial ('WPOT' and 'WPOTNOKICK'; Opp enheimer & Dave1r2006h 
or the circular velocity, v c = VGM n > /j? vl > of the halo from whic h 
the wind is launched (' WVCIRC ; lOppenheimer & Pavel r2008l) . 
In 'WVCIRC the wind velocity and mass loading are given by 
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v w = (3 + n)v c l V2 and 77 = -4= x (v c /v ait ) ', where n and v clit 
are parameters, set to 2 and 150 km s _1 , respectively. In 'WPOT- 
NOKICK' , the wind velocity is given by v w = 3cr, where a is 
the velocity dispersion, calculated from the gravitational poten- 
tial: cr = V-®/2, and 77 = 150 kms"'/(r. In 'WP07" an extra 
velocity kick of 2 X cr is added to each event. These models are 
not directly comparable to those of I O ppenheimer & Pavel I2006L 
l2008l l2009) because these authors used different par ameter values 
and d ifferent models in different papers (see §4.9 of ISchave et al.l 
l201Ch . Moreover, contrary to Oppenheimer & Dave, we do not de- 
co uple the wind particles from the hy drodynamics. As wa s shown 
bv lDalla Vecchia & Schavel d2008h and lSchave etatl ( l2010h . and as 
we will show in Sec. 14. 51 temporarily turning off pressure forces for 
wind particles makes the winds more efficient and has a large effect 
on galaxy properties. 

One important difference between the models discussed here 
and those in Sec. 14.21 is that here the total amount of energy put 
into the wind per unit stellar mass formed is not constant, but de- 
pends on the mass of the halo in which it takes place (E oc in 
'WVCIRC'), and we caution that the total energy used in feedback 
exceeds the total available energy from SNe for the most massive 
galaxies^. Most worryingly, the amount of momentum in the winds 
also exceeds the momentum available from SNe and radiation, as- 
suming that every photon is absorbed once to drive the outflow (so 
no boost fro m optical thickness in the IR), by a factor of about 7 
jHaasll2010l see also the discussion in iDave et alj|20lTh . It should 
therefore be kept in mind that these models may inject an unphysi- 
cally large amount of momentum. Inspection of Fig.[2]reveals that 
these wind prescriptions are capable of almost completely disrupt- 
ing the galaxy disk. 

We compare the effects of the different momentum-driven 
wind models in Fig. [6] All the simulations studied in this section 
predict relatively flat relations between M tot and / star (panel A), 
yielding a flatter galaxy mass function at the low-mass end (panel 
I). None of the models show a characteristic halo mass above which 
feedback becomes ineffective (panels B - F), indicating that the 
feedback remains effective at all masses. In all of these simulations, 
the combination of a large mass-loading in low-mass haloes and a 
high wind-velocity in high-mass haloes means that feedback from 
star formation is always able to drive strong winds and suppress SF. 

4.5 Hydrodynamically decoupled winds 

Many SPH simulations that use kinetic SN feedback employ a 
method known as decoupling that allows SN-driven winds to ef- 
fectively escape galaxies before they begin to interact with the 
halo gas. In these simulations, wind particles, once launched, are 
temporarily decoupled from the hydrodynamics until they escape 
the ISM. During decoupling a gas particle experiences gravity, but 
feels no hydrodynamic drag. Decoupling prevents shock-heating 
and hence the radiative losse s which may otherwise quen ch the 
outflow in high-pressure gas toalla Vecchia & Schavell2012h . and 
it prevents the wind from entraining any surrounding ISM. For a 
detailed study of the effect of decoupling for the case of i solated 
disk galaxy simulations, see iDalla Vecchia & Schavel d2008h . 

To explore the effect of decoupling, one of the OWLS simula- 
ti ons, 'WHYDRODEC , employ s decoupling, following the recipe 
of lSpringel & Hemquistl {2003a). Every particle that is kicked into 

2 This occurs for M tot > 10 12 5 M Q , although the exact mass of equality is 
redshift dependent, due to the redshift dependence of the virial radius. 



the wind feels no hydrodynamic forces for either a fixed amount of 
time (50 Myr), or until the density of the wind particle falls below 
some value (10% of the star formation density threshold, i.e. when 
»h < 10~ 2 cm -3 ). If the wind would retain its original velocity, this 
would correspond to a travelling distance of roughly 30 kpc, so it 
ensures that SN-driven winds escape the galaxy. 

Inspection of Fig. [2] reveals that decoupling the SN-driven 
winds dramatically increases the density of gas around this galaxy 
as every wind particle is able to escape from the galaxy (see also 
IDalla Vecchia & Schavdr200lil) . Decoupling the winds means that 
SN feedback remains capable of suppressing SF to significantly 
higher halo masses, until the wind velocity falls below th e escape 
speed from the halo (see also Springel & Hernquist 2003b). This is 
visible in Fig. [7] where the galaxies in the l REE' simulation show 
sharp upturns in M stM (panel A), SFR (panel B), f lSM (panel D) 
at a halo mass ~ 10 1125 M Q , which are not present in the 'WHY- 
DRODEC simulation, indicating that the SN feedback remains ca- 
pable of quenching SF if winds are artificially decoupled from the 
hydrodynamics. The decoupled winds are capable of ejecting more 
gas entirely from both the galaxy (panels D and F) and the halo 
(panels C and E). 

Interestingly, for the lowest halo masses the SFR is higher for 
' WHYDRODEC than for 'REF ? (panel A). This is expected, be- 
cause wind particles will drag other particles into the wind, thereby 
increasing the effective mass loading, provided the initial wind ve- 
locity is sufficiently high for the winds to leave the galaxy. 

The differences between 'WHYDRODEC and 'REF' are very 
significant, up to 0.5 dex in the ISM fraction (panel D). From this 
we can conclude that gravity (which acts on the winds in both sim- 
ulations) is not the process that makes the winds ineffective. It is 
rather the hydrodynamic drag and the associated radiative losses 
that make the winds less able to escape from high-mass haloes (see 
Dalla Vecchia & Schave 2008, 2oT3: ICreasev et alj|201ll) . 



4.6 Thermal SN feedback 

All the models that we considered up till now, implemented feed- 
back from star formation kinetically. In this section we inves- 
tigate what happens if, instead of launching the wind by in- 
jecting kinetic energy, we inject thermal energy into the gas 
surrounding each newly form ed star particle. As described in 
Dall a Vecchia & Schavel J2012h . the 'thermal feedback' is imple- 
mented stochastically. If the SN energy is distributed amongst all of 
a star particle's neighbours, then the rise in the temperature of each 
particle is so low that cooling times remain very short and the parti- 
cle immediately re-radiates all of the energy. In this case, the feed- 
back wil l have little effec t , unless the cooling is temporarily sup- 
pressed jMori et alj|l997l: iThacker & Couchmar]|200Cl: [Kay et al. 
| 2002j : feommer- Larsen et alJl2003l : lBrook et al.ll2004l : IStinson et a! 
2006). Therefore, we choose to inject the thermal energy into 
neighbouring gas particles stochastically by specifying a temper- 
ature jump, AT = 10 7 5 K, and then calculating for each neighbour- 
ing gas particle the probability that it is heated such that the expec- 
tation value for the total injected energy agrees with the amount of 
feedback energy that is available. This method has some similarit y 
with the promotion feedback model of Scanna pieco et al.l d2006t) . 
We do not turn off radiative cooling at any time. The simulation that 
employs this model is termed 'WTHERMAL' . Given our choices 
for Ar and IMF, the expectation value for the number of heated 
particles per s tar particle is about 1.34 for a fully ionised plasma 
(equation 8 of lDalla Vecchia & Schave1l2012l) . We inject ~ 40% of 
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Figure 7. As Fig. [T] but showing only the 'REF' (black, solid curve) and the simulation in which the wind particles are temporarily decoupled from the 
hydrodynamics ('WHYDRODEC ; red, dotted curve). Despite using the same wind parameters in both of these simulations, artificially decoupling the SN- 
driven winds from the galaxies makes the winds very efficient at suppressing star formation. 



the available SN energy in order to facilitate comparison with the 
other models. 

'WTHERMAL' is compared to 'REF' in Fig. [8] Although the 
two models inject the same amount of SN energy per unit stel- 
lar mass formed, 'WTHERMAL' suppresses star-formation less ef- 
fectively than 'REF' at most masses. This is visible in panels 
A-D, where the thermal feedback simulation lies systematically 
above the kinetic feedback model. The differences are small at high 
masses, because there both types of feedback are effective. How- 
ever, at the lowest halo masses (M tot < 10 10 - 5 Mq, close to our reso- 
lution limit), ' WTHERMAL' predicts lower stellar masses, but still 



slightly higher star formation rates, than 'REF' . This suggests that 
the thermal feedback is more efficient than the kinetic feedback in 
the poorly resolved galaxies that are not plotted here. 



Our finding that the thermal feedback is less efficient t han ki - 
netic feedback is consistent with Dalla Vecchia & Schavd J2012h . 
who predict that for Ar = 10 7 5 K, radiative losses should become 
significant for densities nrr <; 1 crcT 3 (see their equation 18), which 
are certainly reached in all but the lowest-mass galaxies in our sim- 
ulations. 
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Figure 8. As Fig.[T] but showing only the simulation in which SN feedback is implemented thermally ('WTHERM'; red, dotted curve) and 'REF' (black, solid 
curve). Both thermal and kinetic feedback suppress star formation, but for the models presented here, thermal feedback is somewhat weaker, except for the 
lowest, and in these simulations poorly resolved, halo masses. 



4.7 AGN feedback 

AGN feedback is implemented using the method of 
iBooth & Schavel (120091) which is, in turn, a modification of 
that of ISpringel et alj d2005l) . We frequently run a halo finder 
and insert a seed mass black hole (BH) (m sce d = 10 5 /r'M Q ) 
into every halo with mass > 4 x 10 10 h~ 1 M that does not yet 
contain a BH. These seed black holes then grow both through 
merging and gas accretion (which is limited to the Eddington 
rate). Accretion rates in low-density gas («h < 10 -1 cm~ 3 ) are 
assumed to be equal to the Bondi-Hoyle rate. For higher-density, 



star-forming gas the Bondi-Hoyle rate is boosted by a factor 
(;ih/10 _1 cm 3 ) 2 to compensate for th e lack of a cold, inters tellar 
gas phase and the finite resolution (see Booth & Schaye 2009, for a 
full discussion). The BH growth rate is related to its accretion rate 
by otbh = (1 - e r )™ a ccr, where e, = 0.1 is the radiative efficiency of 
the BH. The amount of energy coupled to the surroundings of the 



BH is then given by E = eie t m x 



where c is the speed of light 



and ef is a free parameter, the 'feedback efficiency', that is tuned 
to reproduce the global BH density at z = 0. In the fiducial runs 
ef = 0.15. This model reproduces the observed black hole scaling 
relations ( IBooth & Schavel 120091) . and the observed X-ray and 
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Figure 9. As Fig.[T] but comparing only the simulation with AGN feedback ('AGN'; red, dotted curve) to 'REF' (black solid line). AGN feedback effectively 
ejects gas from galaxies in all haloes where black hole seeds have been placed (M tot > 4 X 10 /i~'M , thus strongly suppressing the star formation rate. 



optica l properties of low-redshift galaxy groups dMcCarthv et al.l 
2010). Note, however, that our simulations have a substantially 
higher resolution than the I00h~' Mpc, 512 3 runs used to compare 
with low-redshift observations. As shown by Booth & Schaye 
(2009), the model is not fully converged for intermediate masses 
with higher resolution, resulting in more efficient BH growth until 
the growth becomes self-regulating. This model does not include 
'radio-mode' feedback (which may be necessary to get a sharp 
exponential drop-off in stellar mass function) and we have not 
varied the parameters of this model in this paper. 

Fig.|2]shows that AGN feedback is the most destructive form 
of feedback, and in a massive galaxy at z = 2 the AGN destroys the 



gaseous disk. This effect is also visible in panels C - E of Fig. [9] 
where it is apparent that the ISM and halo gas content (and thus the 
baryonic content, which is dominated by gas in the halo) of galaxies 
is strongly suppressed relative to the 'REF' simulation. 

In the very lowest-mass resolved haloes (M tot < 10' 06 M Q ), 
'AGN' and 'REF' show virtually identical results because in these 
haloes BHs have not yet been seeded. However, even slightly above 
these halo masses, the BH accretion prescription used in these sim- 
ulations is very efficient at this resolution and the BHs grow and ef- 
fectively suppress star formation (panels A, B, G and H of Fig.[9). It 
has been argued that AGN feedback is necessary to suppress rapid 
cooling of hot halo gas and suppress star formation in high-mass 
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haloes (e.g.lDi Matteo et alJl2005tlCroton et al.ll2004lBower et al.1 
120081 : Booth & Schavdl2009l : iMcCarthv et alJl2010h . In panel (B) 
of Fig. [9] we show the halo SFR as a function of mass for 'AGN' . 
The effect of AGN is indeed very strong at high masses. Although 
the galaxy sSFRs are almost an order of magnitude below the ob- 
servations (panel G) this may not be a problem. The observations 
are of galaxies that are actively star forming. While in the 'REF' 
simulation all galaxies are active, in the 'AGN' simulation there 
are both active and passive galaxies at all masses, leading to a me- 
dian sSFR that lies below the relation obtained for active galaxies 
alone. The stellar mass function (panel I), slightly undershoots the 
observed stellar mass function in the range where both observations 
and simulations exist. This may occur because the SN feedback in 
the 'REF' simulation was tuned to reproduce the peak in the global 
star formation rate density of the Universe. Including AGN feed- 
back will then under-produce the amount of stars. 



5 CONCLUSIONS 

We have analyzed a large set of h igh-resolution cosm ological sim- 
ulations from the OWES project (Sc have et alj|201(il) . focusing on 
the baryonic properties of the galaxy population at redshift 2. In 
particular, we studied the effects of variations in the input physics 
on the stellar mass, SFR, baryon fraction, ISM fraction, and gas 
fraction as a function of halo mass, on the sSFR, ISM mass and the 
gas consumption time scale as a function of stellar mass, and on 
the galaxy stellar mass function. In this paper we focused on varia- 
tions of the parameters of the sub-grid models for radiative cooling, 
feedback from star formation and AGN, as well as the box size and 
the resolution. In Paper II we concentrate on the remaining sub- 
grid prescriptions (e.g. those for star formation), which tend to be 
less important for galaxy properties other than the mass and density 
distribution of the ISM. 

A central conclusion from this work is that the SFR in star 
forming galaxies is self-regulated by galactic winds driven by mas- 
sive stars. Many of our results can be understood if the SFR adjusts 
so that the (time averaged) rate at which energy and momentum are 
injected is sufficient to drive a large -scale outflow that, on average, 
balances the gas accretion rate (e.g. lSchave et alJ|20ld : lDave et al.l 
l2012h . Chemical feedback is also important, because, for a fixed 
redshift and halo mass, the accretion rate is sensitive to the radia- 
tive cooling rate and hence to the metallicity. 

Our results for variations of the kinetic feedback parameters 
and radiative cooling rates can be summarised as follows: 

• Feedback from star formation is only effective if the initial 
wind velocity used for kinetic feedback is sufficiently high. The 
required wind velocity increases with the mass of the galaxy and 
also if metal-line cooling is included. If wind particles are tem- 
porarily decoupled from the hydrodynamics, then the feedback re- 
mains effective up to much higher masses. All these results can 
be understood if radiative l osses inside the ISM are resp onsible for 
the quenching of the winds jDalla Vecchia & Schavehoil. see also 
ICreasevetalj201lh . 

• If the winds are not quenched, then the SFR is inversely pro- 
portional to the amount of energy that is injected per unit stellar 
mass. For a fixed initial wind velocity, this rate of energy injection 
is determined by the initial wind mass loading. 

• If, on the other hand, the initial wind velocity is too low for 
the winds to escape the galaxies, then the behavior of the simu- 
lations tends to that of a simulation without any feedback at all: 
catastrophic cooling resulting in excessive star formation. 



The main conclusions drawn from other variations are: 

• Using a top-heavy IMF for star formation at high pressure 
mainly influences the simulated halos through the extra amount of 
SN energy available per unit stellar mass. 

• Feedback can remain highly efficient at all masses if the initial 
wind velocity is increased with the galaxy mass while keeping the 
momentum that is injected per unit stellar mass constant, as mo- 
tivated by models of winds driven by radiation pressure (although 
the momentum-driven winds used here and in the literature may 
use much more momentum than is actually available). 

• AGN feedback is very efficient at reducing the star formation 
rate and gas content, especially at high masses. 

We compared our predictions to three different observational 
results: the molecular gas mass as a function of stellar mass, the 
sSFR as a function of stellar mass, and the galaxy stellar mass func- 
tion. The latter can be thought of as a convolution of the halo mass 
function and the stellar mass as a function of halo mass. The com- 
parison with observations revealed that: 

• For most simulations, the mass in the ISM in simulated galax- 
ies is slightly lower than the observed molecular gas masses as a 
function of stellar mass. Note though, that the observe d gas masses 
may have been overestimated dNaravanan et alj[2012l) . 

• Except for models with inefficient feedback, the galaxy stellar 
mass function is close to the observed one over much of the ob- 
served mass range. The shape is different though, with most simu- 
lations predicting a steeper low-mass end than (Schechter-like ex- 
trapolations of) observational results. Models with higher initial 
wind mass loading factors predict shallower faint-en d slopes, as 
appea rs to be required by the observations (see also iBower et al.l 
I2OI2I) . None of the simulations predict a clear exponential cut-off 
at the high-mass end, but this could be due to our limited box size 
or lack of 'radio-mode' AGN feedback. 

• The predicted sSFRs as a function of stellar mass tend to be 
lower than observed. For the observed mass range the discrepancy 
is worst if the feedback is efficient. The high observed sSFR can 
only be matched if the feedback is inefficient at the observed mass, 
but highly efficient at lower masses. The observed negative slope in 
the relation between the sSFR and stellar mass is only reproduced 
for galaxies for which feedback is inefficient. 

Thus, there is tension between the observed galaxy stellar 
mass function and the observed sSFRs. The high observed sS- 
FRs are difficult to match unless feedback suddenly becomes in- 
efficient at the stellar masses for which observations are available 
(M, > 10 9 " 5 Mo). It can, however, not remain inefficient as the stel- 
lar mass increases or else the galaxy mass function would become 
too high at the high-mass end. 

We have shown that winds driven by feedback from star for- 
mation determine the main properties of galaxies residing in haloes 
of a given mass (the scatter among the red lines in Fig.[T]is much 
larger than for the blue lines that represent the simulations appear- 
ing in Paper II). Even for a fixed amount of feedback energy per 
unit stellar mass, variations in the sub-grid implementation, e.g. 
different wind velocities and mass loadings (that can be either con- 
stant or fixed functions of local physical conditions), provide us 
with considerable freedom to 'tune' galaxy properties. This free- 
dom can be exploited to match observations spanning a wide range 
of masses, which would provide the simulations with some of the 
attractions of semi-analytic models. However, this potential success 
comes also with the disadvantages of such models: the underlying 
physics may remain poorly understood. As higher resolution simu- 
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lations become feasible, the freedom provided by subgrid models to 
generate galactic outflows in cosmological simulations will be re- 
duced, as the results b ecome less sensitive to the ma nner in which 
the energy is injected jPalla Vecchia & Schavell2012h . 

Further improvement in our understanding of the physics that 
determines the global properties of galaxies will likely come from 
theoretical models and observations focusing on galactic winds. 
The physics of star formation is less crucial as self-regulation im- 
plies that the time-averaged, galaxy-wide SFRs are determined by 
the large-scale inflows and the efficiency with which star formation 
drives galactic winds. 
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Figure Al. Median SFR as a function of halo mass at z = 2 for 5 simu- 
lations with different particle numbers and/or box sizes as indicated in the 
legend. The vertical dotted lines indicate the mass of 2000 dark matter par- 
ticles in the simulations shown by the curves in the corresponding colours. 
At the low-mass end, the median SFR falls to zero, as more than half of the 
haloes in a bin do not have gas particles with a density above the star forma- 
tion threshold. Above a mass corresponding to 2000 dark matter particles 
per halo, the SFR as a function of halo mass is reasonably well resolved. 
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Figure A2. The median sSFRs of haloes as a function of their stellar mass 
at z = 2 for 5 simulations with different particle numbers and/or box sizes 
as indicated in the legend. The vertical dotted lines indicate the mass corre- 
sponding to 100 star particles in the simulations shown by the curves in the 
corresponding colours. The sharp cut-off at low masses again stems from 
the fact that there is a minimum to the (non-zero) SFR. Right of the vertical 
dotted lines the sSFRs are reasonably well converged. 



APPENDIX A: NUMERICAL CONVERGENCE 

Here we explore the sensitivity of our results to numerical resolu- 
tion and simulation box size, and assess down to what mass limits 
we can consider our results numerically converged. All simulations 
used in this section employ the same physical model as the 'REF' 
simulation, but are run with different particle numbers and box 
sizes. We denote the simulations with ' LXXXNYYY' , where XXX 
is the size of the simulation box in co-moving /j~'Mpc and YYY is 
the number of particles per spatial dimension (for both dark mat- 
ter and baryons we use YYY 3 particles). In this nomenclature, the 
'REF' simulation is 'L025N512' . In order to independently exam- 
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Figure A3. Median stellar mass fraction as a function of halo mass at z = 2 
for 5 simulations with different particle numbers and/or box sizes as indi- 
cated in the legend. The vertical dotted lines indicate the mass correspond- 
ing to 2000 dark matter particles in the simulations shown by the curves 
in the corresponding colours. The diagonal black dotted line indicates the 
relation for haloes with 100 star particles (the cut that is made in the rest 
of the paper where relations with stellar mass are shown). As can be seen, 
in the highest resolution simulation, the cuts made throughout this paper in 
minimum number of dark matter particles and the minimum number of star 
particles roughly correspond to each other. At lower resolutions, the cut in 
dark matter particle number is more stringent. 

ine the effects of box size and resolution, we compare two sets of 
simulations: 

• L025N512, L012N256 and L006N128: These simulations have 
the same resolution, but the box size is varied in steps of a factor 
of two from 25 Mpc/A to 6.25 Mpc//7. These three runs are shown 
as black curves with different line styles (solid, dashed and dot- 
dashed, respectively) in all figures in this appendix 

• L025N512, L025N256 and L025N128: These simulations have 
the same box size, but different resolutions. The dark matter particle 
mass in these simulations is 6.3 x 10 6 M G //i, 5.1 x 10 7 M B /h and 
4.1 x 10 8 M Q /h, respectively and the maximum proper gravitational 
softening is 0.5 kpc/h, 1 kpc//j and 2 kpc//i, respectively. In all 
figures in this appendix, these three simulations are shown as solid 
black, blue and red curves, respectively. 

All plots in this appendix include all haloes identified by the 
FoF algorithm (the lowest mass haloes contain 20 dark matter par- 
ticles). Fig. lAll shows the median SFR as a function of halo mass. 
It is immediately clear that simulation box size has no influence 
on the SFRs, as the lines with different line styles (L025N512, 
L012N256 and L006N128) show almost perfect agreement. The 
only effect of box size is that a larger box allows one to sample 
more massive, rarer objects. 

The simulations are less well converged with respect to reso- 
lution (compare the solid black, red and blue curves in Fig. lAlb . 
The vertical, dotted lines denote 2000 times the dark matter par- 
ticle mass in the simulations of the same colour. The halo SFRs 
are reasonably converged above these halo masses, so we employ 
2000 dark matter particles as our resolution limit when comparing 
galaxy properties as a function of halo mass. Observe that the halo 
mass regime where the median SFR is zero is effectively removed 
when demanding a minimum number of 2000 dark matter particles 
per halo, because more than half of the haloes with a low number 



of particles do not have any gas particles with densities above the 
star formation threshold. 

The build-up of stellar mass is influenced by the SFR at all 
epochs prior to the epoch at which it is measured. As all haloes were 
initially small and thus poorly resolved, the early build up of stellar 
mass is underestimated. We therefore expect the convergence of the 
(s)SFR as a function of stellar mass to be worse than that of the SFR 
as a function of total halo mass. 

Fig. IA2I shows the sSFR as a function of stellar mass. The 
vertical cut-off at the low-mass end corresponds again to haloes 
for which the median SFR is zero. The mass range over which the 
sSFR appears approximately converged with respect to resolution 
corresponds to ~ 100 star particles (indicated by the vertical dotted 
lines). Convergence can be seen by comparing the solid black and 
blue lines rightwards of the the blue, dotted line and by comparing 
the solid blue and red curves rightwards of the red, dotted line. We 
note that, as expected, the same trends are found for SFR/M sas (not 
shown). 

Finally, in Fig. [A3] we show the median stellar mass fraction 
as a function of halo mass. The vertical dotted lines indicate our 
adopted resolution limit of 2000 dark matter particles. The diag- 
onal dotted lines indicate the stellar mass fraction for haloes con- 
sisting of 100 star particles, which is our resolution limit for plots 
with stellar mass on the horizontal axis. The fact that for a given 
resolution (i.e. line colour), the solid curve intersects the two dot- 
ted lines in about the same place, implies that the cuts of 100 star 
particles and 2000 dark matter particles are comparable for these 
simulations. Above the limit of 2000 dark matter particles, the stel- 
lar mass fractions are nearly converged. At lower resolution, the cut 
in dark matter particle number is more stringent than a the cut on 
the minimum number of star particles. Therefore, throughout the 
paper we use a minimum of 2000 dark matter particles when we 
plot quantities as a function of halo mass. 

In summary, our resolution requirements are as follows: To 
avoid biasing the results due to resolution effects, we impose a cut 
of 2000 dark matter particles when looking at relations with total 
halo mass and of 100 stars particles when investigating correlations 
with stellar mass. 



APPENDIX B: THE EFFECT OF HALO DEFINITION 

Throughout this p aper all of the ha lo masses that we quote are 
Friends-of-Friends foavis et al.ll985h masses, meaning that all par- 
ticles associated with a given FoF halo contribute to its mass. This 
particular definition of mass is somewhat arbitrary, so we check 
in this appendix whether using a different halo mass definition 
would significantly affect our results. I n particular, we co mpare FoF 
masses to those returned by SubFind toolag et al.ll2009l) . for both 
the main halos and the subhaloes. SubFind iteratively calculates, for 
each FoF halo, the mass of both baryonic and collisionless matter 
that is gravitationally bound to the same structure. 

In Fig. IB II we show the galaxy stellar mass function in FoF 
(solid curve) and SubFind (dotted curve) haloes. The two stellar 
mass functions are virtually identical. This is largely because the 
stars live preferentially in the centres of haloes, and changing the 
halo mass definition affects mainly particles at the edges of the halo. 
This plot suggests that all of our results that are plotted as a function 
of stellar mass are robust to the definition of halo used. 

The halo mass is somewhat more affected by the mass defini- 
tion. In Fig. lB2l we show the median star formation rate as a func- 
tion of halo mass. The solid curve shows the FoF haloes, and the 
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Figure Bl. Galaxy stellar mass function for the 'REF' simulation with two 
different halo mass definitions. The solid (dotted) curve shows the stellar 
mass function when FoF (SubFind) haloes are used. Results for mass func- 
tions are independent of the halo mass definition used, suggesting that any 
results plotted as a function of stellar mass are robust to changes in the halo 
mass definition. 



Figure B2. Relation between halo mass and galaxy star formation rate for 
the 'REF' simulation with two different halo mass definitions. The solid 
(dotted) curve shows the relation for FoF (SubFind) haloes. Although at the 
highest masses, the different mass definitions give results that differ by up 
to 0. 1 dex, this difference is far smaller than the magnitudes of the effects 
we are probing, and also affects all of the different simulations equally. 



dotted curve shows the SubFind haloes. Here we see almost perfect 
agreement at the low-mass end, but the overall effect of SubFind 
at the high-mass end is to unbind haloes that are artificially linked 
together by the FoF algorithm, so at the highest masses the me- 
dian halo mass can decrease by up to 0. 1 dex. We note, however, 
that this same effect will exist in all of the simulations analysed 
here, so although these curves may shift by a small amount, the 
main focus of this paper is on the relative differences between the 
simulations, and these remain totally unaffected by the halo mass 
definition. Note also that the Subfind algorithm only removes mass 
from FoF haloes and never adds any. Many studies use the mass of 
a sphere containing an overdensity of ~ 200 as total mass for the 
main subhaloes, which results in halo masses closer to the FoF halo 
masses used in this paper. 
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